#include "Molmodel.h"
#include <iostream>
#include <fstream>
#include <exception>

using namespace SimTK;
const char* PdbFileName = "SampleRNA.pdb";

int main() {
try{
	//molecule-specialized simbody System
	CompoundSystem system;
	SimbodyMatterSubsystem matter(system);
	GeneralForceSubsystem forces(system);
	DecorationSubsystem decorations(system);

	//molecular force field
	DuMMForceFieldSubsystem forceField(system);
	forceField.loadAmber99Parameters();

	forceField.setCoulombGlobalScaleFactor(0);
	forceField.setGbsaGlobalScaleFactor(0);

	RNA rna("AUG");
	rna.assignBiotypes();
	system.adoptCompound(rna);

	NoseHooverThermostat thermo(forces, matter, 293.15);

	system.addEventReporter(
		new Visualizer::Reporter(system, 0.020));

	std::ofstream pdbFile; pdbFile.open(PdbFileName);
	system.addEventReporter(
		new PeriodicPdbWriter(system, pdbFile, 0.1000));

	system.modelCompounds();

	State state = system.realizeTopology();

	LocalEnergyMinimizer::minimizeEnergy(system, state, 15.0);

	VerletIntegrator integ(system);
	TimeStepper ts(system, integ);
	ts.initialize(state);
	ts.stepTo(10.0);

	pdbFile.close();
	std::cout << "\nWrote pdb file "
              << Pathname::getCurrentWorkingDirectory()
			  << PdbFileName << "\n";
    return 0;
} 
catch(const std::exception& e) {
    std::cerr << "ERROR: " << e.what() << std::endl;
    return 1;
}
catch(...) {
    std::cerr << "ERROR: An unknown exception was raised"
              << std::endl;
    return 1;
}
}